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Abstract 

We investigate the continuum q-Potts model at its transition point from the dis- 
ordered to the ordered regime, with particular emphasis on the coexistence of 
disordered and ordered phases in the high-g case. We argue that occurrence of 
phase transition can be seen as percolation in the related random cluster represen- 
tation, similarly to the lattice Potts model, and investigate the typical structure 
of clusters for high q. We also report on numerical simulations in two dimensions 
using a continuum version of the Swendsen-Wang algorithm, compare the results 
with earlier simulations which used the invaded cluster algorithm, and discuss 
implications on the geometry of clusters in the disordered and ordered phases. 

KEYWORDS: Continuum Potts model, first-order phase transition, continuum 
percolation, continuum Swendsen-Wang cluster algorithm 



1 Introduction 



The Potts model is one of the classical models of statistical mechanics exhibiting 
a phase transition. In its standard version, it is denned on the square lattice Z 2 , 
where it was first studied ^B]- The parameter controlling its phase transition is 
temperature: at sufficiently high temperatures the Gibbs measure is unique, while for 
low temperatures the number of translation invariant extremal Gibbs measures (pure 
phases) coincides with the cardinality of the spin state space (usually denoted by q) 
|Sj. Moreover, when q is large enough, these regimes meet at a specific value of the 
temperature where q distinct "ordered" phases and one "disordered" phase coexist. 
This was shown in i!3 j by using reflection positivity arguments, and in [21 by 
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Pirogov-Sinai theory. The transition is in this case of first order meaning that exactly 
at the transition point the derivative of the pressure with respect to temperature has 
a jump. 

The mechanism of phase transition is well understood. When q is in the low range, 
the appearance of phase transition is explained by ground state degeneracy, similarly 
to the Ising model: the Ising model has two ground states (indeed, it is equivalent 
with the Potts model for q = 2), while the Potts model has exactly q such states. As 
q is taken larger, there is ever increasing possibility for entropy to take domination, 
and for large enough q there is a temperature at which energy and entropy balance 
each other, leading to coexistence of the high and low-temperature phases. 

To understand the structure of these pure phases, one has to have a notion of what 
a typical spin configuration looks like in each such Gibbs state. Presently, there is a 
suggestive picture offered by the random cluster representation. According to this, 
phase transition in the Potts model occurs exactly at that value of the temperature 
for which the associated random cluster model is at the percolation threshold 
llUj . This is possible to prove for d > 2, and in fact for quite general underlying 
lattices. The structure of the percolation clusters allows a good insight into the 
typical configurations of the Potts pure phases. 

In this paper we consider the similar problem of non-uniqueness for the Potts 
model placed in the continuum instead of a lattice. In this case, the particles are not 
sitting in lattice points but at random points in IR d . The a priori distribution of the 
set of occupied positions is Poisson process with intensity z, which is then modified 
by the Potts interaction — a repulsion between particles of different types. (In the 
case of zero temperature, this model coincides with the multitype hard-core exclusion 
model of Widom and Rowlinson US]-) The crucial parameter is then the activity z: 
It is known that, for each number q of different types and at any fixed temperature, 
there is only one Gibbs measure whenever z is small enough, while for high enough 
values there are q translation invariant extreme Gibbs measures; see [21 El ESI an d the 
more recent papers (31 El using random cluster methods. However, in the continuum 
setting it has not been proved so far that these two regimes meet at a single critical 
activity z c , and that the phase transition at z c is of first order when q is large enough. 

The difficulty lies in the fact that there are no obvious extensions to the continuum 
of the methods available for the lattice Potts model. On the one hand, reflection 
positivity does not apply to continuum models. On the other, the contour techniques 
of Pirogov-Sinai theory used in [2"lll4| do not appear to admit an immediate extension. 
These papers either deal with the continuum Widom- Rowlinson model, but then only 
cover the case when z is large and q is held fixed (not necessarily large), proving 
coexistence of just q phases, or show the coexistence of q+1 phases at z c for large q, 
but this only for the lattice Potts model. 

As a rigorous proof is still lacking, it should be worthwhile to report on further 
progress. Based on the random-cluster representation of the model, we investigate the 
structure of clusters in terms of their dissociation probabilities under resampling of 
locations, and show that clusters with positive dissociation probabilities are unlikely 
to occur uniformly in the activity z when q is large. (We believe that these ideas will 
eventually lead to a rigorous proof of first-order phase transition, but at this stage 
there are still a number of difficulties to overcome.) On the other hand, we have 
undertaken a numerical study of these problems. We use a natural continuum analog 
of the Swendsen-Wang cluster algorithm which was originally designed for the lattice 
Potts model. Arguing that a coexistence of ordered and disordered phases manifests 



2 



itself as a medium-term dependence of the algorithm on the initial conditions, we find 
that for any q there is only one critical activity z c , and the phase transition at z c is 
of second order when q = 2,3,4, while it is of first order for q > 5. 

This picture confirms the main results of earlier numerical studies |12[ 1161 l2*T] 
based on the so-called invaded cluster algorithm. However, the earlier results were 
not conclusive about the order of the transition in the q = 4 case which we can better 
determine now. In addition, the exact relationship between the stationary measure of 
the invaded cluster algorithm and the corresponding finite volume Gibbs measures of 
the Potts model is not known. It has only been postulated that in the infinite volume 
limit there should be a unique activity z such that they coincide, but this remains 
still to be proven. 

Therefore, we felt it necessary to make an independent numerical study using an 
algorithm whose stationary measure we can show to be a Gibbs measure of the Potts 
model. This enables us to check, at least numerically, some of the results obtained by 
the invaded cluster algorithm. As already mentioned, our results confirm qualitatively 
the earlier ones, and it is plausible that the invaded cluster algorithm describes the 
continuum Potts model in the infinite volume limit. However, we also found a case in 
which the invaded cluster algorithm gave results with significant finite size bias, slowly 
decreasing with the volume. Such behavior makes the control of finite size effects in 
the invaded cluster algorithm difficult, and a careful analysis of the bias would be 
advisable before the results are actually applied to describe the Gibbs states of the 
Potts model. 

The outline of this paper is as follows: In Section 2 we introduce the continuum 
Potts model and the continuum Swendsen-Wang algorithm. In Section 3 we analyze 
the typical structure of clusters in the high-g continuum Potts model in terms of their 
dissociation probabilities, and discuss the behavior of the algorithm in the presence of 
first-order phase transitions. Section 4 provides our simulation results, while Section 5 
contains a detailed discussion including a comparison with the results obtained by the 
invaded cluster algorithm, and a discussion on the structure of pure phases. 

2 Model and algorithm 

2.1 Potts model in the continuum 

The continuum Potts model is a model of point particles having q > 2 different 
types and sitting in a rectangular box A C M. d , d > 2. Rather than of particles 
of different types, one may also think of particles with a ferromagnetic spin with q 
possible orientations. A configuration of particles in A is given by a pair X = (X,a), 
where X is the set of occupied positions, and a : X — > {l,...,g} is a mapping 
attaching to each particle in X its type, or "color". Writing X a = {x 6 X : a(x) = a} 
for the configuration of particles of type a, we may also think of X as the g-tuple 
of the pairwise disjoint sets X a belonging to Xa = {X C A : #A < oo}, the set 
of all finite subsets of A. The configuration space is thus equal to X A q \ the set of 
g-tuples of pairwise disjoint elements of Xa- The particles are supposed to interact 
via a repulsive interspecies pair potential <p : M d — > [0, oo] of bounded support. For 
simplicity we confine ourselves to the case of a step potential 
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already considered, e.g., in |15[ 120] . The Hamiltonian in A is thus given by 

f a (x)= Yl E <p( x -v)- (2- 1 ) 

l<a<fe<<j x£X a ,y£X b 

Here we impose periodic boundary conditions, meaning that the difference x—y has 
to be understood modulo A. (We note in passing that one could also add a molecular, 
type-independent interaction term, as was done in [7j- Here, however, we stick to 
the simple case above.) The associated Gibbs distribution with activity z > and 
temperature T > is then 

<? 

/iAAT(dX) = ZxXt e M~ h aW/T] J! z#Xa LA(dX a ) , (2.2) 

o=l 

where La is the Lebesgue-Poisson measure on X\ defined by 

L A (A) = V )— / 1a({xi, . . . ,x n }) dxi---dx n 

n! /in 

n>0 JA 

and ^a.,z,t is the normalizing constant. The zero-temperature case corresponds to the 
classical model of Widom and Rowlinson (23] with hard-core interspecies repulsion. 
One can then imagine that the particles form balls of diameter 1 which can overlap 
only when they are of the same type. 

The behavior of the infinite-volume Gibbs states of this model is completely un- 
derstood when z is either small or large. If z is small enough, there is a unique 
infinite-volume state which is disordered, in that it is invariant under permutations 
of particle types; this can be seen, for example, by using disagreement percolation; cf. 
Proposition 7 of 0. If z is sufficiently large (depending on T), there exist q distinct 
phases which are ordered, or demixed, in that one particle type is more frequent than 
all other types; see (3J E3 for the case T = and El fo r general T. It is expected, 
but not rigorously known, that there exists a sharp activity threshold z c = z c (T) such 
that the infinite- volume Gibbs state is unique when z < z c and non-unique for z > z c . 
(This lack of knowledge is due to the fact that the model does not have any useful 
stochastic monotonicity properties. The only monotonicity known is that the particle 
density is an increasing function of z; cf. Section 4.2 of [Sj and equation (|4.1|) .) If q is 
large enough, it is further expected that the transition at z c is of first order, meaning 
that the disordered and the q ordered phases exist simultaneously. This is the problem 
we address in this paper. 

2.2 Random-cluster representation 

Just as the lattice Potts model, the continuum Potts model admits a random-cluster 
representation of Fortuin-Kasteleyn type; see jH) and the references therein, as well as 
OIZ]. This random-cluster representation will become important in the following. The 
random-cluster measure associated to (|2.2j) is a probability measure for random graphs 
r = (X, E) in A. The vertex set X is obtained from the configuration X = (X, a) 
by disregarding the particle types described by a, and the edge set E is obtained by 
drawing random edges between the points of X. Specifically, for each X £ X A let 
£x consist of all sets of non-oriented edges between pairs of distinct points of X, and 
vx,T be the probability measure on Ex for which an edge between a pair {x, y} C X is 



4 



drawn, independently of all other edges, with probability p(x—y) = 1 — e~ v ^ x ~ y ^ T ; as 
before, the difference x—y is understood modulo A. (In the Widom-Rowlinson case 
of hard-core interspecies repulsion, the randomness of the edges disappears in that all 
points of distance < 1 are connected automatically.) 

The random-cluster measure associated to (12. 2 1) thus lives on the space Q\ = \T = 
(X, E) : X £ X\, E £ £x} of all finite graphs in A, and is given by 

XA,z,T(dX, dE) = Z K \ T z * x q k{x ^ L A (dX) v XyT {dE) , (2.3) 

where k(X, E) stands for the number of clusters of the graph (X, E), and ^a,z,t again 
denotes the normalization constant. (Note that this definition makes sense for any 
real q > 0.) As indicated by our notation, the normalization constant is in fact the 
same in either of equalities Q2.2JI and l|2.3|) for any allowed values of the parameters. 
This was established in |J| as part of the proof for the following precise relationship 
between the two measures: 

Proposition 2.1 (p rx x) Take a particle configuration X = (X,a) £ xjf 1 ' with 
distribution ha,z,t and define a random graph (X, E) £ Q A as follows: Independently 
for each pair {x,y} of points of the same type (i.e., a(x) = cr(y)) let {x,y} £ E with 
probability p(x—y) = 1 — e~' fi ( x -y)/ T , Then (X, E) has distribution Xk,z,T- 

(X ^ A*) Pick a random graph T = (X, E) £ Q\ according to XA,z,T and define 
a type assignment a as follows: For each cluster C of T assign a type a £ {l,...,q} 
independently and with equal probability, and then define cr(x) = a for all x in the 
union of all clusters of the type a. Then X = (X, a) has distribution fJ,\ Z T- 

To obtain a joint picture of the continuum Potts model and its random-cluster 
representation, one should think of cluster-colored graphs T a = (X,E,a), where 
r = (X, E) £ Q\ and a is a mapping attaching to each cluster C of T a color 
cr(C) £ {l,...,q}. Let us use the notation (/)a,z,t for expectation values of a 
random variable / on the cluster-colored graphs. The continuum Potts model is 
then obtained by interpreting a as a function on X which is constant on all clus- 
ters, and then forgetting the edges E; that is, for any / which depends only on 
X = (X, a), (/)a,z,t = / /M.,z,T(dX) /(X). Likewise, forgetting the colors one arrives 
at the random-cluster measure: for / depending only on (X,E), we have (/)a,z,t = 
f X A,z, T (dX,dE)f(X,E). 

2.3 Conditional single-type distributions 

The Gibbs distribution ha,z,t in <|2.2j) has also another useful property easily to be 
exploited for simulation. Namely, if we fix all particles except those of a given type 
a £ {1, . . . , q}, then the conditional distribution of the particles of type a is Poisson 
with a simple intensity function. Specifically, for any bounded non-negative function 
u on A let 

nl(dX) = exp [ - / u(x) dx] JJ u(x) L A (dX) (2.4) 

A xex 

be the Poisson point process on X\ with intensity function u. Then the following 
observation follows immediately from the definitions (|2,lj) and (|2.2I) : 

Proposition 2.2 Let 1 < a < q be a given type and, for any X = (X,a) £ 

let X^ a = {x £ X : a(x) ^ a} be the set of positions of all particles having types 
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different from a. Then, under /m.,z,t> the conditional distribution of X a given X+ a 

is equal to the Poisson point process w ith intensity function zp(- \X^ a ), 

where 

p(x\X^ a ) = exp - (p(x-y)/T . 

To simulate tvV' X ^ a one can use the well-known fact that X ^ a can be 

obtained from the homogeneous Poisson point process 7r^ with constant intensity 
function z by a random thinning: each point x from a 7r^-sample X is kept, indepen- 
dently of all other points, with probability p(x\X^ a ); otherwise x is removed. In the 
spirit of the random-cluster representation discussed above, this can also be achieved 
by independently drawing (virtual) edges between the points x of X and y of X^ a 
with probability p(x—y), and deleting all x G X that are connected by an edge to 
some y G X^ a . 

2.4 The continuum Swendsen-Wang algorithm 

The algorithm of Swendsen and Wang [22] is by now a standard device for simulating 
the lattice Ising and Potts models. It can be characterized as the algorithm which 
alternatively applies the transition probabilities relating the Potts model with its 
random-cluster representation. The naive analog for the continuum Potts model would 
be an alternative application of the two steps described in Proposition 12.11 Note, 
however, that these steps always keep the set of occupied positions fixed. That is, 
these transition steps are unable to equilibrate the particle positions. (Iterating these 
steps, one would rather arrive at the discrete random-cluster distribution of edges 
between the vertices chosen initially.) So, one has to combine these steps with a 
further simulation step which takes care of the positions. The simplest such step is 
the Gibbs sampler based on the conditional probabilities of Proposition 12.21 We are 
thus led to the following continuum version of the Swendsen-Wang algorithm, variants 
of which have already been proposed independently in [Ij and : 

Continuum Swendsen-Wang algorithm: Start from any initial configuration 
X G and iterate the sweep consisting of the following three steps: 

CSW 1: Resampling of positions. Successively for a = 1, . . . , q, replace X a by a sample 
from the Poisson point process 71^ \ x ¥=a) ^ us i n g a random thinning of 7T^. 

CSW 2: Drawing edges. Let X = (Ja=i X a and, independently for each pair {x, y} of 
points of the same type, draw an edge from x to y with probability p(x—y) = 
1 — e~^ x ~ y ^ T . Let E be the resulting set of edges, and consider the graph 
T = (X,E). 

CSW 3: Choice of types. For each cluster C of T, independently of all other clusters, 
pick a random type uniformly in {1, . . . , q} and assign this type to each x G C. 
Let X a be the set of vertices receiving type a, and X = (X\, . . . , X q ). 

In view of Propositions 12.11 and 12.21 it is clear that the Gibbs distribution ^a.z,t is 
invariant under this algorithm. In fact, the following ergodic theorem holds: 

Proposition 2.3 Let X n G n > 0, be the realization of the continuum Swendsen- 
Wang algorithm after n sweeps. Then the distribution of X„ converges to ^k,z,t in 
total variation norm at a geometric rate. 
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Proof: It suffices to observe that, in Step CSW 1, X a = for all a with probability at 
least 5 = e _2g ' A L So, if X n and X^ are two versions of the process starting from dif- 
ferent initial configurations but otherwise using the same realizations of randomness, 
then 

Prob(X n + X' n ) < (1 - 5) n , 
whence the proposition follows immediately. □ 

For practical purposes, particularly in our context, the ergodic theorem 12. 31 is rather 
misleading. This is because the rate of convergence towards ^a,z,t can be extremely 
small, even for A of moderate size. (This is already seen from the number 5 above, 
though this is only a simple lower estimate of the coupling probability.) In particular, 
this is the case in the presence of a first-order phase transition when fi\zT is essentially 
supported on disjoint sets Ai, i = 0, . . . , q, that are typical for the q + 1 coexisting 
phases. In this case, the sets A dls = Aq and A ord = Uf =1 Aj are separated by tight 
bottlenecks of the CSW-algorithm. In fact, over a fairly long initial period the CSW- 
algorithm will converge to the conditional probability /iA z t(* I A), with A = A dls 
or A ord depending on the initial condition, and [aa,z,t is reached only after a time 
that is far exceeding any reasonable observation period. In this way, one can detect 
a first order phase transition by comparing the CSW-algorithm for different initial 
conditions. We will discuss this point in more detail in Subsection 13. 21 

We conclude this section comparing the algorithm described above with the related 
algorithm invented in [1J . 

Remark 2.4 Instead of the systematic scan through all types in Step CSW 1 above, 
one could also use a random scan by resampling only the positions of a random (or, 
by type symmetry, a fixed) type. Contracting our Step CSW 3 with the successive 
Step CWS 1 (for a single a) one then arrives at the following algorithm proposed in 

lanu: 

Random-scan continuum Swendsen-Wang algorithm: Starting from any initial 
graph T £ Qa, iterate the following two steps: 

1. Each cluster C of T, independently of all others, is deleted with probability 1/q 
and retained with probability (q — l)/q. Let (X \d, S Ld) be the remaining graph. 

2. Choose a sample X new from 7r A P ^ ' Xold ^ and, independently for each pair {x, y} 
in X ncw , draw an edge from x to y with probability p(x—y). Let E new be the 
resulting set of edges, and consider the graph T = (X \d U X ncw , E \d U E new ). 

While this version has the advantage of remaining completely in the random-cluster 
picture (and thus working also for non-integer q), a priori it is not evident whether it 
is more efficient or not. Clearly, q sweeps of the random-scan version require the same 
numerical effort as one sweep of the systematic-scan version. However, the former 
has higher correlations even after q sweeps since the waiting time until all clusters 
are resampled (the maximum of the geometric waiting times for replacement of a 
single cluster) has expectation larger than q, and thus does not seem to converge at 
a geometric rate not depending on the number of clusters. We performed a brief 
numerical comparison of the efficiencies and behavior of the above two algorithms, 
and these support the above picture: both algorithms lead to similar behavior but 
the systematic scan is slightly more efficient than the random scan. 
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3 Detecting first-order phase transitions in finite volume 



3.1 The structure of clusters for high q 

In this subsection we present a rigorous result showing a common feature of all clusters 
at any activity when q becomes large, in arbitrary dimensions. We recall first what 
happens in the planar lattice Potts model; see |13j . For large q it is known that, with 
probability one, the plaquettes (2x2 squares) on which the configuration shows one 
particular of several typical patterns form an infinite cluster. In the ordered phase with 
dominating spin value a the pattern corresponds to the local ground state in which 
spins take the same value a, while in the disordered phase the pattern corresponds to 
a local ceiling state in which all nearest neighbor spins differ. Configurations which 
belong to neither category fail to get weight in the thermodynamic limit. In the 
continuum Potts system we expect a similar characterization in terms of percolation: 
an ordered phase with dominating type a should be characterized by percolation of 
spins of value a, while the disordered phase by percolation of vacancies. Moreover, 
it should be possible to derive these properties from certain typical local patterns 
characteristic of the ordered or disordered case. Corollary 13.21 below will show which 
kind of local patterns can occur for any z and any phase when q is large. 

In the following we confine ourselves to the Widom-Rowlinson case T = 0. (We 
believe that similar estimates should hold also for T > 0, but this would need extra 
effort.) Consider a configuration X G X\ in a box A. Since we are in the Widom- 
Rowlinson case, the associated set of edges in the random cluster model is determin- 
istic, viz. Ex = {{x,y} C X : \x — y\ < l}. Let C C X be a cluster of the graph 



T = (X,E X ). We write 

U(X \ C) = {x G A : 3 y G X \ C, \y - x\ < 1} 

for the part of A in which any point is connected to X\C, and Aq(X) = A\U(X\C) 
for the available free space of C. We consider the probability 



is the distribution of a configuration of N particles thrown independently and uni- 
formly into A (where |A| is the Lebesgue measure of A), and £;(£) is the number of 
clusters of (£, £^). We call 5c(X) dissociation probability, for it measures how big 
the chance is to split C into two disconnected parts by a random resampling of its 
points in the room available after taking C away. In particular, a small dissociation 
probability makes it unlikely for C to admit a pivotal point which cannot be removed 
without splitting C into disconnected parts, and therefore expresses some kind of ro- 
bustness of C. The following result states that, for large q, all clusters are robust in 
this sense, and this is a property uniform in z. 

Proposition 3.1 Let T = and k(A) = m&x{k(X) : X G X\} be the maximal 
number of clusters in A. Then for all q > 0, 5 > and z > we have 



S C (X) = L Ac{x)]#c (C : k(0 > 2). 



(3.1) 



In the above, #C is the number of particles of C, 
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Note that /c(A) is of the order |A| when A is a square. Therefore, this estimate cannot 
be applied directly to the infinite volume limit. The main point of the bound is that 
it is uniform in the activity z, and thus applies to all possible finite- volume "phases" 
for sufficiently high values of q. We have not tried to optimize the constant in the 
Proposition, and 2 K ( A ) is most likely very far from being optimal. 

Postponing the proof for a moment, let us first dicuss the consequences of this 
result. Intuitively, a weak dissociation tendency of a cluster C means that either C is 
a singleton (in which case it cannot dissociate), or the available free space Ac(X) is 
small compared to the number of particles. One instance of the latter case is captured 
by the following definition: For any given < 7 < 1 , let a cluster C in a configuration 
X be called ^-confined if #C > 2 and Ac(X) admits no two Borel subsets Ai and A2 
such that they are at least a unit distance apart and each contains a fraction 7/2 of the 
total free volume, i.e., such that dist(Ai,A2) > 1 and |Ai|,|A2| > 7I Ac(X)|/2. Of 
course, this condition means that Ac(X) must be small, in that it either has diameter 
at most one or consists of a solid core exceeding the diameter 1 only by some tiny 
filaments (possibly scattered all through A). 

Corollary 3.2 Let T = 0, < 7 < 1, a large number Nq S N, and e > be given. 
Then there exists some qo > 1 such that, for all q > qo and z > 0, 

XA A o(4 is U<)>l-e, 

where A° A ld = {X e X A : 3 cluster C C X, #C > N } and 

Af s = {x eX A :V clusters C C X, #C = 1 or C is 7-confined} . 

Evidently, the event A A describes a scenario characteristic for the disordered 
phase: For very small z, all clusters will typically be singletons, and the model is 
similar to the hard-core gas of balls of unit diameter. If z increases, the singletons 
(considered as balls) will become more and more densely packed. At the close-packing 
density of balls the color entropy cannot be increased any more, and the increasing 
particle density will force the system to build up confined clusters of two or more 
overlapping balls. But all clusters are still separated by channels of vacancies. At the 
threshold z c , color entropy breaks down in favor of positional entropy, which means 
that the system will form a large cluster C of size #C > Nq, where particles have more 
positional "degrees of freedom". It is then plausible that Aq(X) fills a macroscopic 
part of A and all other clusters are confined. Since the color is constant on C, we 
could conclude that a fixed spin value a is dominating with overwhelming probability. 
So, with a proper choice of A^, the A° A rd scenario should be typical for the ordered 
phases above z c . 

We now turn to the proofs of Proposition 13.11 and Corollary 13.21 

Proof of Proposition We start by noting a symmetry property of the Lebesgue- 
Poisson measure L\ which readily follows from its definition: For any measurable 
function F : X A — > [0, oof, the expression 

I L A (dX) £ f L A | #€ (dr?) F(e, V,X\0 (3.2) 
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is invariant under the exchange of the first two arguments of F. Now, 



Sq Z Kzfi XKzfi(x G X A : 3 cluster C C X, S C (X) > 6^ 

< f L A (dX) Z# X q k ^ +1 1« cluster of X} 6((X) 

< / L A (dX) l{fcA € (jr),fc(o=i} 

x (|A|/|A 5 (X)|)#€ / L^Wl^^^aj^^)^^) 



since + 1 < k(X \ £ U rj) under the circumstances described by the indicator 

functions. Next we use the symmetry property of expressions of the form (j.3.2j) . 
together with the fact that Ag(X) depends only on X \ £. The last integral then 
becomes 

/ L A (dX) £ / 

L A\M( dr l) 1 { V cA ! .(X),k( v )=l} 

J £CX J 

x (|A|/| A f (X)|)#" l { ^ Adx) , m > 2} z* X q k ^ 

= f L A (dX)z# x q k ^ Y, htcA^HOm L Adx)m {k(-) = 1) . 

Finally, we estimate away the last probability in the last integrand simply by 1 and 
note that the condition £ C Ag(X) means that £ is disconnected from X \ £ and, 
therefore, consists of a union of clusters of X. Since there are at most 2 k ( x ^ < 2 K ( A ) 
such unions of clusters, the last expression is not larger than Z a<z< t 2 k ^ , and the 
result follows. □ 

Proof of Corollary El Let 5 = j N °/2 and q be so large that 2 K ^/5q < e. By 
Proposition 13 . 1| it will be sufficient to show that 5c (X) > 5 whenever 2 < j^C < Nq 
and C is not 7-confined. However, in this case there exist two Borel subsets Ai and 
A2 such that dist(Ai, A2) > 1 and |Aj| > 7|Ac*(X)|/2 for i = 1,2. So, a resampling 
of the points of C within Ac(X) will certainly produce at least two clusters whenever 
all new points fall within Ai U A2, and each of the sets Ai and A2 gets at least one 
of them. Hence 



1 



>S, 



and the proof is complete. □ 



3.2 Quasi-absorbing sets of the CSW-algorithm 

We now take up the discussion begun after Proposition 12.31 on the behavior of the 
CSW-algorithm in the presence of first-order phase transitions. We will argue that 
a first-order phase transition is characterized by the appearance of two different sets 
which are nearly absorbing for the CSW-algorithm, so that its behavior over a rea- 
sonable observation period strongly depends on the initial condition. The occurrence 
of such a dependence on the initial condition is therefore an indication of a first order 
phase transition. 
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Our reasoning consists of three parts: First we will ask how a first order transi- 
tion will manifest itself in finite volume, then we will study the influence of a quasi- 
absorbing set on the medium-term behavior of the CSW-algorithm, and finally we 
explain why a jump of the particle density should imply a bottleneck in the CSW- 
algorithm. 

1. How can a first-order phase transition be observed in finite volume? By def- 
inition, a first order transition is identified by a discontinuity of the first derivative 
of the infinite volume pressure. In our case, with z as parameter, this corresponds 
to a discontinuity of the particle density at the transition point z c . It is conjectured 
that for the continuum Potts model with sufficiently large q, the percolation threshold 
z c should drive such a first order transition. In addition, at z c one expects a coex- 
istence of disordered and ordered phases, in that there exist two mutually singular 
type- invariant Gibbs measures, |U dls y anc ^ Z^^T' wn i cn can be obtained as limits of 
the unique type- invariant measures for z f z c resp. z I z c , and are the only extremal 
elements of the set of type-invariant Gibbs measures. In particular, the average par- 
ticle density of /x dls T should be different from that of fM^ d T - Since the measure (aa,z c ,t 
is also type-invariant, its infinite volume limit would then be a (presumably non- 
trivial) convex combination of these measures. There are then two disjoint sets of 
configurations A^ 1S and A™ d , approaching the disjoint supports of HzT an d A^'r m 
the infinite volume limit in the sense that ha,z c ,t{A^ 1s U A^ rd ) — > 1 as A f M. d , while 
separately each set has a probability strictly bounded away from zero, and the con- 
ditional measures (J,\,z c ,T 

(• | 4 is/ord ) have average densities which stay a fixed value 
apart from each other. Since in both cases the variance of the density tends to zero in 
the infinite volume limit, we can also assume the sets to be chosen so that the density 
distributions are are almost mutually singular. 

2. What is the behavior of a Markov chain admitting a unique invariant measure fi 
almost concentrated on two disjoint sets A\ and A2 such that the conditional measures 
iu(- 1 Ai) can be distinguished by an observable /? It is then plausible to expect that 
the sets Ai are nearly absorbing, in that the Markov chain stays within these sets 
with probability close to 1. The following remark explores the medium-term behavior 
of / for such a Markov chain. 

Remark 3.3 Let (X n ) ra >o be a Markov chain with a Polish state space E and sta- 
tionary distribution /u, A a measurable subset of E with fJ,(A) > 0, ha = h('\A) 
the associated conditional distribution and (X^) n >o the induced Markov chain on A 
with invariant distribution ha, and suppose (X^) n >o is uniformly ergodic. Also, let 
/ : E — > M rf be any measurable observable (not constant on A) and e > 0. By the 
large deviation principle for Markov chains there exists then some 5 > such that 



/ii N r 

Prob bE/( x «)- / f d ^ 

^ n=l J 



>e\< e~ m 



for all N > 1; see Theorems 6.3.8 and 2.3.6 of jSJ. Finally, suppose that A is nearly 
absorbing for (X n ) n >o, in that 

Prob(X n+ i G A I X n ) > 7 when X n G A , 

where 7 is so close to 1 that < for some prescribed error probability a. 

It then follows that, for large but not too large N, the time average of / along the first 
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N steps of the original Markov chain (X n ) ra >o is close to the conditional expectation 
J f d/iA, with probability close to 1. Namely, 

Prob(|-i^/(X n )- ! fdfi A \>e) < 2a when X e A and < AT < . 

^ n=l ^ ' 

Indeed, (X n ) n >o coincides with (X^) n >o up to the first time T such that (X n ) n >o 
leaves A, and T is geometrically distributed with parameter 1 — 7. Splitting the event 
in question into two parts according to whether T < N or not we therefore obtain the 
upper bound 1 — 7^ + e~ SN , and the result follows. 

In other words, if a MCMC algorithm has a very small probability of leaving a 
set A, it typically stays for the whole observation period within A and thus, during 
this period, coincides with the induced Markov chain on A which converges to fj,A- In 
the case under consideration, it turns out that the mixing properties of the induced 
Markov chains are good enough for this scheme to work for both a set A = A^ d of 
ordered configurations, and a set A = A^ s of disordered configurations. We cannot 
yet give a proof for the existence of such absorbing sets for the values of q we inspected 
here, but we believe they have a characterization similar to what was given for the 
high-g limit in Corollary 13.21 (whence we have chosen the same notation). In fact, the 
simulations described below support the existence of such absorbing sets and make 
clear that a first-order transition of the continuum Potts model does indeed manifest 
itself by a pronounced "bottleneck" of the CSW-algorithm between two different quasi- 
absorbing sets. 

3. Finally we ask: Are there any indications that the CSW-algorithm does indeed 
have a bottleneck between two different quasi-absorbing sets? We will argue that, for 
any given e > 0, a jump greater than e in the particle density in one sweep of the CSW- 
algorithm has probability tending to zero in the infinite volume limit (apart from the 
very first sweeps). Therefore, if the sets ^4^ rd and ^4^ ls are defined in a way allowing a 
distinction by different particle densities then, near a first order transition point and 
for large enough volumes, these sets should be nearly absorbing, and passages from 
one to the other should create a bottleneck after a large enough number of sweeps, as 
this would require to pass through a gap in the density distribution of ^k,z c ,t- 

The density of a configuration can only be changed in step CSW 1. Suppose 
that the input configuration is (X, a) £ X-jfi . We remove all particles of type a, and 
consider the available free volume, denoted by Ax a (X) as in Sec. 13.11 For the Widom- 
Rowlinson model, the new particles of type a are distributed exactly according to a 
Poisson measure with activity z in this volume. Since the density distribution given 
by the Poisson measure on a set with Lebesgue measure v has expectation value z and 
variance z/v, the new set has a density z with fluctuations of the order of \jyfv. If 
v -C |A|, these changes to the total particle density are negligible, and if v 3> 1, then 
the change in the total particle density is sharply concentrated to (zv — #X a )/\k\. 
When T > 0, particles can also be added to the complement of Ax a (X), but even 
then the probability is significant only where the old particles are not too dense, that 
is, typically only near the boundary of the set Ax a (X). 

If this is not the first sweep, then X a was obtained by the same procedure and 
should consist of regions which either are small containing not too many particles, or 
have particle density close to z; see also the discussion in Sec. 13. II for the high-g limit. 
Then the above argument indicates that typically the particle density should fluctuate 
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like a Poisson density, that is C(|A| -1 / 2 ), in one sweep. Our numerical results, to be 
discussed in the next section, agree with this scaling relation, apart from second order 
phase transition points. 

4 Numerical results 

4.1 Implementation of the algorithm 

In this section we present our results from applying the CSW algorithm to the two- 
dimensional Potts model. The main numerical hurdle to overcome in the simulation of 
the algorithm is the large number of particles in the kind of volumes we need in order 
to get reliable estimates about the properties of the infinite volume phases. Here the 
bottleneck is not the memory required to store the configurations, but rather the time 
required to find those "old" particles which lie within the interaction radius from a 
given "new" particle, as well as the time needed for dividing the particles into clusters 
according to a given edge configuration. 

The second part, the forming of clusters for a given configuration of particles and 
edges, can be done very efficiently by using the algorithm described in ^7]. We used 
the tree-based union/find algorithm given in section II B. of JJj, while simultaneously 
keeping track of the size and of the "corners" of the clusters (see section 14.21 for the 
precise definition). 

To overcome the large volume problem in the first part of the algorithm, we used 
"hashing" of the box into smaller cells, altogether N? of them. The number of cells 
in one dimension, was chosen so that the average number of particles in the 
cell, as determined by the Poisson process with activity z, would be about 10. For 
each cell, we created q directed lists, one for each possible type. The list number a 
contained pointers to all those particles which had the type a and which were within 
the interaction radius (here the unit distance) from the cell. 

This information was used to substantially reduce the time needed both in creation 
of the thinned Poisson configuration, and in computation of the open bonds: by 
construction, if we add a particle of type a anywhere in the cell, then it can interact 
only with particles in one of the lists of the cell with a' ^ a. The use of q separate 
lists allowed for easy removal of the particles with a certain type which was needed 
in the first part of the algorithm. 

All simulations were performed using square boxes of linear size L, i.e., A = [0, L] 2 , 
with periodic boundary conditions. We employed three different initial conditions: a 
Poisson sample with activity z, either assigning all particles the type 1, or choosing 
the types randomly, and a "disordered crystal" where the particles lie in a certain 
dense square lattice with alternating types. The initial conditions with a uniformly 
colored Poisson sample are called here "ordered". The other two alternatives - ran- 
domly colored Poisson and the disordered crystal - led to the same behavior (i.e., to 
measurements within error bars of each other) in all those cases where we tried both. 
Therefore, we call them collectively "disordered" initial conditions in the following. 

4.2 Measurements 

For measurement of the properties of the infinite volume Gibbs states, we employed 
the numerical CSW-update algorithm with several values of the box size L to an initial 
configuration Xq £ ■ After a preset number of steps, no > 1 (called "burn- in" or 
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equilibration period), we measured the value of an observable /(r CT ) for the following 
consecutive jim > steps, and obtained a sequence of samples f n = f(T° +n ) for 
n = 1, . . . ,riM- The values no and um were assumed to be chosen so large that the 
average would well approximate the corresponding expectation value, 



-I "M 

Efn ~ (/}A,z,T- 



nM 

n=l 

Finding a good choice for no and nM was not straightforward. It was particularly 
difficult near the phase transition points where we found out that, even for this cluster 
algorithm, the equilibration and decorrelation times for relevant observables can be 
very large. We chose, quite arbitrarily, no = 250 and nM = 2500 as a first guess, and 
increased these values when necessary; when quoting the results we will use the short- 
hand phrase "using no + nM sweeps" to give the actual values used in computation of 
the results. 

To estimate the error arising from the finiteness of nM, we computed the standard 
deviation from new samples obtained by dividing the data into 10 blocks: the block 
averages are less correlated than two consecutive samples, and as long as nM > 10 x 
(decorrelation time) this should yield a fairly reliable estimate of the error. Explicitly, 
the "errors" given later were obtained by defining, for «b = ^m/10 and k = 1, . . . , 10, 

^ riB, 

fk = — ^2 f (fc-l)n B +n> 
n=l 

then computing the sample variance 

I 10 l 10 2 1 10 10 

S ) = 9 [12^) 2 - Y^(Y1 /*) ] = 9 - To E /. 



k=l k=l k=l 



of (fk), and estimating the standard deviation of the average of (f n ) by Sf/ylO. Some 
of the decorrelation times for large boxes were indeed very long (see, for instance, 
Figure IHJ) and these more elaborate methods were required to get a sensible error 
estimate. 

For any given colored cluster configuration T a = (X, E, a) we considered the 
following four observables p, p', 7 and dp erc : 

1. Particle density 

p{T°) = N(X)/L\ 

where N(X) = j^X is the total number of particles. 

2. Slope estimator p' . Since 

3 11 L 2 

-p-(p)K,z,T = ~{Np) Kz , T - -(N) A ^ T (p) A ^ T = — Var(p), (4.1) 

c z z z z 

the scaled sample variance of p, 

L z 1 i~ nM 1 / nM 

p = t E p% — [J2p< 

- T — 1 Lf— ' nM V, 



z „ M -n- n M fc=1 



estimates the derivative of {p)\,z,T with respect to z. This is not an observable in the 
previous sense, so we computed its error estimate by using as samples the 10 sample 
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variances computed from the 10 sample blocks of p. Note also that if p' remains 
bounded when L — > oo, then the corresponding standard deviation of density in one 
CSW-step is 0(1/L) near stationarity. 

3. Largest cluster size 7. This observable measures the ratio of particles in the 
largest cluster: 

7 (n = max(#C)/iV(X), 
o 

where the maximum is taken over all clusters C of V = (X, E). 

4- Percolation radius d perc . This quantity measures the spread-out of the clusters 
from a given L-independent set So in the middle of A. As reference set we used 

S = { x e A : 3y such that \x - y\ < 1/2 and \\y - (L/2,L/2)||oo < 3/2} , 

that is, a central 4x4 square with rounded corners. For i = 1, 2, let &r(r) and bf(T) 
denote the minimum and maximum of the coordinates of clusters with particles in So, 
i.e., 

6,~ (r) = min minxj, bj (T) = max maxxj, 
C-.CnSo^i x&c c-.CnSo^Q xeC 

where C runs through the set of all clusters of T. Our definition for the percolation 
radius then reads 

d pcrc (n = max{| - b-(T) - 1, 6+(r) - \ - l,o}. 

(This particular choice is adapted to the Widom-Rowlinson case, T = 0, where the 
most natural percolating objects are the discs of radius 1/2, as explained in Sec. 12.21 
For convenience, we retained this definition also in the case of positive temperatures 
where it can appear to be unnecessarily complicated.) For all practical purposes, it is 
safe to think of d pe rc as the maximal distance the clusters of V percolate away from 
the central 4x4 square. 

After determining the critical values of z, we also repeated some of the simula- 
tions near these values in order to find out how the particles are distributed between 
clusters of different size. To this end, we built histograms for cluster sizes by using 
the observables 

where C(x) denotes that cluster of T which contains the particle x, and the second 
sum goes over all clusters of T. They describe the portion of the particles in clusters 
with size (relative to N) in the interval A. Here the intervals were chosen by dividing 
[0, 1] into 100 pieces, i.e., using = [k — l,/c)/100, for k = 1, . . . ,99, and A100 = 
[99/100,1]. We also measured the average ratio of particles in very small clusters, 
with sizes from 1 to 100. 



4.3 Computation of the critical activity 

For the computation of the critical activity for a fixed temperature T we started from 
a box with side length L = 8 and computed the above observables for several values 
of z by using an equidistant grid in a suitable range, always for both disordered and 
ordered initial conditions. This allowed an inspection of the effect of initial conditions 
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Figure 1: The measured values of p as a function of z for q = 10 and T = using five 
different box sizes and 250 + 2500 sweeps, all with both ordered (black) and disordered 
(grey) initial conditions. Only the equilibrated results have been shown, see the text 
for explanation. 



and, apart from a neighborhood of a first order transition, these values always agreed 
within the computed error bars. 

It was expected that some percolation property could be an order parameter for 
the transition, and it turned out that both d pe rc and 7 had a pronounced change at 
the transition. Using the values measured for p, p', 7 and d perc we could locate the 
critical value z c approximately. The simulations were then repeated in a neighborhood 
of this value on a finer grid, but with twice the length L. This was repeated until 
sufficient accuracy was achieved, typically at L = 128, although we had to go up to 
box sizes L = 512 for q = 4 and q = 5. Figure ^ shows the results of such an iteration 
for the density in the case q = 10, T = 0. 

The order of the transition was determined from the dependence of the observ- 
ables on using either ordered or disordered initial conditions. If the different initial 
conditions led to different values of density, the transition was determined to be of 
first order. Our results also fully support the discussion made in section 13.21 which 
allows us to identify the two different results as properties of the different coexisting 
phases. Figures Q and the q = 5 part of present typical examples of the behavior in 
these instances: both initial conditions lead to density fluctuations 0(1/L), but the 
average values are separated by a constant 0(1). (This is true only for large enough 
L. For very small L, there is a significant probability to jump from one density re- 
gion to another, and then the result is some average of the values for each phase, 
and the standard deviation is of the order of the gap between these values.) In addi- 
tion, the observed values are right-continuous for the ordered initial conditions, and 
left-continuous for the disordered ones. 

In case the values agreed, we next looked at the behavior of p': since its maximum 
was then always found to diverge as some power of L in a neighborhood of the perco- 
lation threshold, we call these second order transitions. In Figure El we have plotted 
the evolution of the density under our algorithm for the values in the borderline cases: 
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Figure 2: The evolution of p/z under the algorithm for q = 4 (the left figure) and 
q = 5 (the right figure), and for two different initial conditions as described in the 
text: ordered (the black line) and disordered (the grey line). In all of these runs, 
T = 0, L = 512, and z ss z c (z = 2.05125 for q = 4 and z = 2.168 for q = 5). n SW 
denotes the number of "sweeps" performed, and the dotted line represents the chosen 
equilibration cutoff, no- 
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2.051(4) 
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2.273(7) 


2nd 
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2.1675(13) 
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2.424(2) 


1st 


10 


2.56(1) 


1st 


2.965(25) 


1st 


50 


3.65(30) 


1st 


4.95(65) 


1st 



Table 1: Estimates for the critical activity z c and the order of the transition in the 
L — > oo limit for several q and at zero and one non-zero temperature. The error 
estimates are fairly conservative, see the text for how they were obtained from the 
simulations. 



the largest q for which the transition was found to be of second order, q = 4, and the 
smallest q for which it was of first order, q = 5. 

In some cases, especially when using the "wrong" initial conditions near the borders 
of a first order transition region, the equilibration times turned out to be much longer 
than the chosen no - see Figure 01 for a typical instance. Similar problems, combined 
with very long decorrelation times, plagued the q = 4 simulations in large boxes 
as well; note, for instance, that no = 250 is clearly insufficient for the left part in 
Figure El However, since typically only one of the initial conditions suffered from 
these long equilibration times, we decided, instead of redoing the simulations with 
very large no, to throw away the non-equilibrated value and use only the equilibrated 
one. This explains some of the apparently "missing points" in Figure ^ Nevertheless, 
we always kept at least one result for each z, redoing the simulations with larger no 
when necessary. 

The results from this analysis are given in Tabled For estimating finite-size ef- 
fects, i.e., the difference of critical values from the L — > oo limit, we used certain 
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Figure 3: The evolution of p/z for q = 10, T = 0, and z = 2.54, using L = 128 and 
ordered initial conditions, n^w denotes the number of sweeps performed, and the 
dotted line the preset equilibration cutoff. 



monotonicity properties observed while doing the simulations. For second order tran- 
sitions, the observable p' exhibits a divergence near the transition points, and the 
value quoted is the position of the maximum of the peak for the largest box size used. 
As can be seen also in Figure 0J this appears to increase monotonously with L, and 
therefore it can reasonably be taken as a lower bound for the limiting value. The error 
given, 5z c , is a value such that for z > z c + 5z c the measurements of p' for the two 
largest boxes agree with each other. Again, as seen in Figure 0] this is quite robust 
a value, fairly independent of L, and would in these cases be better understood as 
an upper bound for the error. For first order transitions, we used the property that 
the "coexistence window" between the disordered and ordered phase goes to zero as 
L — » oo, apparently monotonously. In these cases, the error is given by the range of 
values in which both initial conditions were equilibrated but yielded differing results, 
plus one grid spacing. For instance, in the table we have z c = 2.56(1) for q = 10, 
T = 0, which was obtained from the single coexistence value shown in Figure ^ 

5 Discussion 

5.1 Comparison with earlier results 

The Widom-Rowlinson and continuum Potts models have already been studied nu- 
merically before in ^2] an d jU- These simulations used the so-called "invaded cluster" 
(IC) algorithm introduced in for studying the lattice Potts models near criticality. 
This algorithm is similar to the Swendsen-Wang type of algorithms presented here. 
The main difference is that instead of generating samples for the finite volume Gibbs 
measure directly, a similar updating algorithm with a suitably chosen "stopping rule" 
is used for generating samples for a measure which differs from any of the finite vol- 
ume Gibbs measures, but which is claimed to approach the correct Gibbs measure in 
the infinite volume limit. 

The main advantage of the IC algorithm is that an advance scanning of the pa- 
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Figure 4: The measured values of p' as the function of z for four different box sizes, 
computed for q = 3 and T = using 250 + 2500 sweeps. The value shown is an 
average of the results with ordered and disordered initial conditions. The shaded 
regions depict the results given in ^2] for the critical z with L = 160 (left region) and 
with L = 40 (right region). 



rameter space for finding the transition values is not necessary, as the stopping rule 
is assumed to be chosen so as to force the system to be at the transition point in the 
infinite volume limit. Two stopping rules are used in j!2| I21j to study the continuum 
models: the percolation rule for q < 3 when the transition was expected to be of 
second order, and the fixed density rule for q > 3 to study whether the transition 
is of the first order. The main disadvantage of the IC algorithm is that the above 
mentioned convergence to the correct limit measure has not been proved so far and, 
in any case, it is quite difficult to control the finite size effects. 

In Figure we have compared the results for the critical activity z c for q = 3 and 
T = from the present algorithm to those obtained using the percolating IC algorithm 
in ^2]. The finite size effects appear to be more prominent in the percolating IC 
algorithm: by Figure 0] increasing the linear size fourfold from L = 40 to L = 160 
does not appear even to halve the systematic error to the L — » oo value which we 
estimated (using lattice sizes up to L = 256) from the position of the peak in p'. 

For q = 2, T = 0, we found for the limiting value z c = 1.718(7), while in ^2] 
z c = 1.7262(4) for L = 160, and z c = 1.7201(7) for L = 40 were measured. The finite 
size effects seem to be weaker in this case. In [2^, simulations were performed also 
for a few non-zero temperatures. Unfortunately, there is only one instance in which 
we can directly compare our results with theirs: for T = 0.5, q = 2 we estimated 
z c (L = oo) = 1.86(4) while in |21j for L = 20 with the same parameter values 
z c = 1.8508(4) is given. Due to coarseness of our result, we cannot really compare the 
finite size effects in this case. 

Let us offer a possible explanation for the observed differences between the results 
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from these two methods. Since percolation (i.e., the existence of a cluster spanning 
the whole of A) was used as the stopping condition in the above quoted results, 
every sample configuration contains a percolating cluster. Even if we assume that 
samples approximating the infinite volume measure at the percolation threshold are 
generated by the IC-algorithm, this sampling method introduces some bias into the 
measurements. Actually, comparing the IC results with ours, it appears that the IC- 
method overestimates the critical activity z c . This can be understood by observing 
that, in finite volume, the probability of having a percolating cluster is a continuous 
function of z ascending from to 1 around z c , and thus is certainly not close to 1 at 
z c but only when z is sufficiently larger than z c . 

For first order transitions, we would expect exactly the opposite to happen: the 
percolating IC-algorithm should underestimate z c . Indeed, since this algorithm pro- 
duces only samples with a percolating cluster, the ordered phase gains an advantage 
over the disordered phase whenever there is a chance for it to occur. But, in finite 
volume, the ordered and disordered phases can coexist throughout a whole range of 
parameter values z (rather than only at z c , as in infinite volume). The threshold 
detected by the IC-algorithm therefore identifies only the lower end of the coexistence 
interval. Unfortunately, we cannot test this hypothesis, as the fixed-density stopping 
rule, and not the percolation one, was used in [2^ for obtaining the critical activity 
for those q exhibiting first order transitions. 

Apart from these differences, our results confirm those found in and |21j . 
For both the Widom-Rowlinson model (T = 0) and the Potts model (at least with 
this particular non-zero temperature) we found only a single phase transition point, 
and the onset of percolation is an order parameter for this transition. As in these 
references, we also found that the transition is of second order for q = 2, 3, 4, and of 
first order for q > 5. 

5.2 Structure of pure phases 

Apart from localizing the critical activity z c and clarifying the nature of the tran- 
sition, our simulation measurements also provide some insight into the structure of 
the pure phases. Let us start by considering Figure |S] which shows the evolution of 
our observables p (density), 7 (largest cluster size) and d peTC (percolation distance) 
for q = 5 during the simulation steps for two initial conditions: ordered (left-hand 
side) and disordered (right-hand side). It is clear that p and 7 display more or less 
stationary fluctuations around some value that depends on the initial condition. This 
indicates the stability of the ordered and disordered phases over the observation pe- 
riod and allows to infer that these phases coexist; recall the discussion in Section 13.21 
It also shows that the ordered phases have a higher particle density than the disor- 
dered phase, which means that the particle density in infinite volume should have a 
jump at z c . Likewise, the proportion of particles in the largest cluster is nearly 1 in 
the ordered case and nearly in the disordered case, from which we conclude that 
the typical configurations of the ordered phases contain a macroscopic cluster, while 
those of the disordered phase do not. A glance at the evolution of d per c reveals that, 
in the ordered case, the macroscopic cluster typically hits the central 4x4 square. On 
the other hand, in the disordered case we see many "spikes", telling that it can quite 
well happen that one can walk a long distance from the central 4x4 square along 
the random graph, but the corresponding clusters are quite "fragile" and filamented, 
surviving only for a short time. All these effects become more pronounced when q 
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Figure 5: The evolution of the density, largest cluster size, and percolation distance, 
p, 7, and d perc , respectively (see the text for precise definitions). Ordered initial 
conditions were used in the left case, disordered in the right, while otherwise in both 
cases q = 5, T = 0, z = 2.168 ~ z c , and L = 512. 
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Figure 6: As in Figure but for the case q = 10, T = 0, z = 2.56, and L = 128. 
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gets larger, as can be seen from Figure El for the case q = 10. 

By way of contrast, Figures and |H1 show the corresponding structure in cases 
q = 2 and q = 4. It is obvious that essentially no difference can be found between 
the ordered and disordered initial conditions, and that the criticality of z c manifests 
itself only by very large fluctuations. We thus conclude that the phase transition is of 
second order. In fact, it also becomes clear that q = 4 is a boundary case: the portion 
of particles in the largest cluster is typically quite large, and so is the percolation 
distance. This gives the hint that the onset of percolation at z c is quite rapid, and 
that the underlying value of z in Figure |H1 is actually slightly above z c . 

Figure El presents the cluster-size distributions in the second-order case q = 4 and 
the first-order case q = 5, again for ordered and disordered initial conditions. For 
9 = 4, there is again almost no influence of the initial conditions, while for q = 5 there 
is a dramatic difference, and the different phases are separated by a range of values of 
cluster size which were not observed in either phase: those corresponding to the case 
when about half of the particles are in the maximal cluster. Finally, Figure ITU1 shows 
that in both phases the portion of particles in small clusters decays like a power-law 
of the cluster size (note the log- log scale in the figure). However, in the ordered phase 
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Figure 7: As in Figure EJ but for the case q = 2, T = 0, z = 1.72, and L = 128. 
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Figure 8: As in Figure 03 but for the case q = 4, T = 0, z = 2.05125, and L = 512. 
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Figure 9: Histograms for the probabilities of finding a particle in a cluster which 
contains a portion r of the particles. On the left, the data were obtained from the 
runs with q = 4 depicted in Figure |H1 with the shaded area giving the histograms for 
the ordered initial condition and the gray line for the disordered. The right figure 
shows the corresponding results for the q = 5 runs given in Figure [SJ 
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Figure 10: The probability of finding a particle in a cluster of size s. Measured from 
the same runs as the q = 5 part of Figure |HJ again black corresponding to the ordered 
initial condition and grey to the disordered one. 

the decay is clearly faster. 
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